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SUMMARY 

A time-and-space accurate and computationally efficient fully three-dimensional unsteady 
te mp erature field analysis computer code has been developed for truly arbitrary configurations. It 
uses Boundary Element Method (BEM) formulation based on an unsteady Green's function 
approach, multi-point Gaussian quadrature spatial integration on each panel, and a highly clustered 
time-step integration. The code accepts either temperatures or heat fluxes as boundary conditions 
that can vary tn time on a point-by-point basis. Comparisons of the BEM numerical results and 
known analytical unsteady results for simple shapes demonstrate very high accuracy and reliability 
of the algorithm. An example of computed three-dimensional temperature and heat flux fields in a 
realistically shaped internally cooled turbine blade is also discussed. 


INTRODUCTION 

For linear boundary- value and initial-value problems it is computationally more efficient to 
use BEM rather than finite differencing or finite element technique. An additional benefit is that 
the BEM requires a relatively coarse grid and that it can be easily extended to non-linear 
conduction problems via Kirchoffs transformation. Probably the most important fact is that the 
BEM is essentially a non-iterative technique thus making the BEM codes more reliable 
[Dulikravich and Hayes 1988; Dargush and Banerjee 1991]. 


THEORY OF 3-D BOUNDARY ELEMENTS 
The governing equation for heat conduction in the absence of heat generators is 


du 

dt 


= aV 2 u 


( 1 ) 


where a is the thermal diffusivity such that a = k/pc, k is the thermal conductivity, p is the mass 
density, and c is the specific heat As the problem is now time-dependent, the weighted residual 
statement must be integrated with relation to time and space. After Integration by parts twice 
[Brebbia and Dominguez 1989], one obtains 
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where to < x < t and the 9u*/9t term was obtained integrating by parts with respect to time. 
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The corresponding fundamental solution for this equation is 


"* ’ [iJa(L)] M exp 

where D is the number of spatial dimensionality, that is, D ^ 3 for three-dimensional problems 
and r is the distance from the point of integration, xj, to the observation point, xi. The 
fu ndam ental solution satisfies the auxiliary equation 
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Using these relations, one 



l l 

c iU t + JJuq* drdr ■ JJqu* dTdr + Juu* dfl 

t»r tor jj 


( 6 ) 


The last term in the above formula represents die contribution of the initial state at t — to 
Since the fundamental solution is time dependent, one does not need to solve this equation with 
an iterative scheme as is usually done in finite elements or finite differences. Instead, this 
equation can be solved for any time step after a known initial state although small time intervals 
are recommended. In the latter case, the potential at each node within the domain needs to be 
evaluated at the end of each time step in order to determine the initial conditions for the next time 
step. Although the primary advantage of boundary elements is lost (for unsteady problems, the 
discretization of the volume is necessary), the matrix is much smaller than those in finite elements 
and the inversion needs only be performed once if time-independent boundary conditions are 
enforced. 

The region, O, and bounding surface, T, are discretized into y volume cells utilizing 
a total of nodes and Ngp surface elements using a total N nodes. Nodal quantities are 
interpolated linearly across each individual cell or surface panel. One also needs to assume a 
variation in time for the functions u and q. If these functions do not vary significantly with 
respect to time, they may be treated as constant over small time intervals and the time integration 
may be performed stepwise. If greater accuracy is required, these functions may be assumed to 
vary linearly such that 

U (Xi,Xj,t,T) = Ifia. u + LLL u m 

At At K ' 


where the subscript 0 indicate the variable at the previous time level and t is the current time level. 
The time step is defined as At = t - t 0 . Once fully discretized, the boundary clement equation 
may be expressed as 


[H] U = [G] Q + P (8) 

for constant time elements. If linearly varying time elements are used, the boundary element 
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equation becomes 

[H] U + [H 0 ] U Q = [G] Q + [G 0 ] Q 0 + P (9) 

The term U = is the vector of potentials and Q = (qj,q 2 iQ 3 «—» ( lN S p) is a 

vector of fluxes where each term qj contains four potential derivatives corresponding to the comer 
vertices of the Ngp quadrilateral surface panels. In addition, each term in the [G] matrix contains 
four terms corresponding to each qj term. The terms in the [H] matrix are the properly summed 
coefficients corresponding to the potential at a specific node. The boundary element algorithm 
developed for use in three-dimensions with the inverse design code was formulated with 
isoparametric, quadrilateral surface panels and eight-point linearly deformed parallelepiped 
volume cells. Integration in space was done numerically with Gaussian quadrature and time 
integration was performed numerically using the trapezoid rule. Since the fundamental solutions 

contain a singularity at the end of each time step (at the time level x = t), the time integration points 
were clustered strongly near the singularity, x = L 

A transformation from the global (^T|,Q coordinate system is necessary to a localized 

coordinate system defined over the surface of the body or to a coordinate system defined 
for the volume cell of integration. The volume integrals are transformed such that 


Juu* dQ = JJJuu* IJI d£d77d£ 
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where HI is the determinant of the Jacobian matrix of the transformation, that is. 


| J] = det 
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The terms of the G and H matrices for constant time elements may be written as 


GO) 


( 11 ) 
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for the jth surface panel. The subscript m indicates the node number corresponding to the kth 
vertex of the jth surface panel and Nq is the number of surface panels connected to the mth node. 
The function f^ is the kth interpolation function corresponding to the value at the kth comer vertex 
of the jth surface panel. Discretization of the surface is identical to that described in section 2.2. 

If linearly varying time elements are used, the terms of the [G] and [H] matrices are similar and 
are formed into pairs of coefficient matrices due to the initial state and the current time level at the 
boundary. 

dr J dr 


(13) 


q * dTj df 


(14) 


H, - |jj*V 


The vector P indicates the contribudon of the inidal state throughout the domain on the ith 
observadon point and is of the form 


Pi = 2fX u ojk dQj 

j-l k-l Qj 


(15) 


where 0^ is the kth trilinear interpolation function of the linearly deformed parallelepiped cell and 
is of the form 

- i(l±S)(l±r)Xl±C) (1© 

The term Uqjjj is the potential on the kth vertex of the jth parallelepiped volume cell at time x = tQ. 
In discretized form, the fundamental solution contains the term r representing the distance from 
the point of integration on the jth volume cell or surface panel to the ith control point The normal 
derivative of the fundamental solution may be determined as 


t) 


where d is defined by the normal distance from the ith point under consideration to the surface 
paneL If the control point is on the surface of integration, then q* = 0 and the diagonal of [H] 
may be computed implicitly by the application of a uniform potential over the whole domain, 
which will give zero normal fluxes at the boundaries such that 
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NUMERICAL RESULTS 

The accuracy and reliability of the BEM formulation governing three-dimensional, 
unsteady heat transfer problems has been verified versus the analytical solution for a finite length 
rod. Total surface of a cylinder of 0.5 m radius and 1.0 m in length was modeled with 216 nodes 

and 108 surface panels as shown in Figure 1. The thermal diffusivity, a, was specified to be 1.0 
no? s' 1 . The initial temperature of the cylinder was uniformly 0 °C and contained no heat sources 
or sinks. Then, suddenly, the boundary conditions on the cylinder were specified as 100 °C on 
the front face, 0 °C on the back face while the outer radial surface was kept adiabatic. The BEM 
algorithm used constant time interpolation, 3-point Gaussian quadrature for the surface and 
volume integration and a linear variation of die temperature and heat flux along the surface panels 
and volume cells. The analytic solution for this test case corresponded to the one-dimensional 
unsteady heat flow in a finite thin rod. Temperatures were obtained at various time levels and at 
several axial locations and are compared to die analytic results shown in Figure 2. As seen in this 
figure, there is a discrepancy between the numerical and analytic solutions averaging about 6 ®C. 

The geometry of the cylinder was then modeled differently by slightly clustering the surface 
panels and volume cells near the hot end Figure 3 shows the results of the BEM in this case. 


470 



Notice that the average error has reduced to about 4 °G The unsteady BEM algorithm was then 
developed with a linear variation within each time step for die temperatures and heat fluxes. This 
test case used the same slightly clustered geometry with identical boundary and initial conditions. 

A single BEM analysis run using linear time interpolation consumed about 15% more CPU time 
than the BEM formulation with constant time interpolation. It resulted in a solution which 
averaged a 3 °C error compared to the analytic solution. Figure 4 shows the computed 
temperatures from the BEM using linear time interpolation against die analytic solution at several 
axial locations. Figure 5 is a plot of the error between the BEM and the analytic solution. 

The accuracy of the BEM algorithm could not be improved further while using the cylinder 
test case. Instead, a different and simpler geometry was developed and the BEM solution with this 
new geometry showed that much of the error in the previous tests was due to the cylindrical 
geometry. That is, since the BEM uses flat quadrilateral surface panels, the exact geometry of the 
cylinder surface could not be modeled properly. Figure 1 shows a discretized cross section of the 
circular cylinder and clearly depicts the inability of a limited number of flat surface panels to 
properly capture the surface curvature of the cylinder. Besides, elements surrounding the cylinder 
axis are nearly triangular in shape. These situations produce surface and volume integrals which 
behave somewhat singularly. The result is that the integrals are not properly integrated and may 
involve ill-conditioned BEM solution matrices. 

The new geometry, a rectangular box of 0.1 m x 0.1 m base cross section and 1.0 m in 
length, was then used instead of the cylinder. The geometry was divided into 10 axial cells of 
equal size. The entire surface of the box was discretized with 44 nodes (four nodes per each 
section) and 44 flat quadrilateral surface panels (four side panels per each section). The boundary 
conditions on the box were specified as 100 °C on the front face, 0 °C on the back face while the 
side surfaces were kept adiabatic. The BEM algorithm used linear time interpolation, 5-point 
Gaussian quadrature and a linear variation of the temperature and heat flux along the surface panels 
and volume cells. Temperatures obtained with the unsteady BEM algorithm for the rectangular 
box were compared to the analytic solution at several axial locations. Figure 6 illustrates that the 
BEM solution for the box was much more accurate than those for the cylinder. Figure 7 shows 
the absolute value of the error in the temperatures computed using the BEM. These results indicate 
that the BEM generates an error of 0.5 °C with the maximum error below 0.9 °C. 

The unsteady BEM algorithm was then modified to incorporate temperature-dependent material 
properties. Although the BEM solution of the linear heat conduction equation is quite fast 
(^uiring less than 10 CPU seconds for 25 time steps on an IBM 3090 for the rectangular box), 
the addition of temperature-dependent material properties greatly increases the computational 
e: j ^Normally the BEM solution matrices need only to be computed once if the time intervals 
and diffusivity are constant Since the diffusivity is now a function of temperature, the BEM 
solution matrices need to be developed at each time interval using temperatures computed at each 
source point in the surface and volume integrands. 

The same rectangular box geometry and boundary conditions were used to test the 
accuracy of the three-dimensional, unsteady BEM algorithm with temperature-dependent material 

properties. The reference thermal conductivity was Xo = 1.0 kcal nr 1 s* 1 K* 1 and it varied 

linearly with temperature as X = Xq+ C T. The temperature variations in time at a single axial 
location were collected for various degrees of nonlinearity given by the parameter C. These results 
are shown m Figure 8 and compare well with published computational results involving finite 
elements [Tanaka, Kikuta and Togoh 1987]. Total CPU time for 25 time steps with the 
temperature-dependent physical properties was approximately 300 seconds on an IBM 3090. 


CONCLUSIONS 

A fully three-dimensional unsteady heat conduction analysis code has been successfully 
developed and tested against known analytical solutions. The code is computationally efficient 
and reliable and can be used on arbitrary configurations. A modification involving temperature- 
dependent thermal diffusivity was also incorporated and shown to produce good results. 
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Figure 1 . Geometry of a cylinder for the verification of the three-dimensional, 
unsteady BEM formulation. 
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gure 2 


Comparison of temperatures between the unsteady BEM solution 
and the analytic solution using constant time elements. 



Figure 3 . Comparison of temperatures between the unsteady BEM solution 
and the analytic solution using a.refined grid and constant time 
elements. 
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Figure 4 . Comparison of temperatures between the unsteady BEM solution 
and the analytic solution using a refined grid and linear time 
elements. 
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Figure 5 . Relative error in temperature between the unsteady BEM solution 
and the analytic solution using a refined grid and linear time 
elements. 
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Figure 6 . Comparison of temperatures between the unsteady BEM solution of a 
reoanguJar box and the analytic solution using linear lime interpolation 
and accurate quadrature integration. 



Figure 7 . Error in temperature between the unsteady BEM solution of a 
rectangular box and the analytic solution. 
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Figure 8 Results of the unsteady BEM algorithm with temperature-dependent 
^ material properties. The figure shows temperatures versus nme at the 
z = 0.2 m axial location for various ranges of nonlinearity. 
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